home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / fft.lha / fft / fft_init.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  2KB  |  97 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *                  Fast Fourier Transform
  6.  *
  7.  *             Initialization of the FFT packet
  8.  *
  9.  *  The file defines the constructor and the destructor for the FFT class,
  10.  *  allocates global arrays and initializes global data
  11.  *
  12.  ************************************************************************
  13.  */
  14.  
  15.  
  16. #include "fft.h"
  17. #pragma implementation "fft.h"
  18.  
  19. #include <math.h>
  20. #include <Complex.h>
  21. #include <std.h>
  22.  
  23.  
  24. /*
  25.  *------------------------------------------------------------------------
  26.  *            Constructor and Destructor
  27.  */
  28.  
  29. FFT::FFT(const int n, const double dr=1)
  30.      : N(n), dr(dr)
  31. {
  32.   assure( n > 3, "At least 4 points must be given for transforms");
  33.   assure( dr > 0, "Grid mesh in the r-space, dr, must be positive");
  34.  
  35.   register int i;
  36.   for(i=1, logN=0; i < N; i *= 2, logN++)
  37.     ;
  38.   assure( i == N, "No. of points has to be the exact power of two");
  39.  
  40.   A = new Complex[N];
  41.   W = new Complex[N];
  42.   index_conversion_table = new short[N];
  43.   A_end = A + N;
  44.   
  45.   fill_in_W();
  46.   fill_in_index_conversion_table();
  47. }
  48.  
  49. FFT::~FFT(void)
  50. {
  51.   assert( N > 0 && A != 0 && W != 0 && A_end == A + N );
  52.   delete A;
  53.   delete W;
  54.   delete index_conversion_table;
  55. }
  56.  
  57. /*
  58.  *------------------------------------------------------------------------
  59.  *        Fill in the index_conversion_table so that
  60.  *  if j = J[m-1] * 2^(m-1) +  J[m-2] * 2^(m-2) + ... + J1 * 2^1 + J0 * 2^0
  61.  *  then index_conversion_table[j] =
  62.  *       J0 * 2^(m-1) + J1 * 2^(m-2) + ... + J[m-2] * 2^1 + J[m-1] * 2^0
  63.  *
  64.  *  m being Log2(N)
  65.  */
  66.  
  67. void FFT::fill_in_index_conversion_table()
  68. {
  69.   register int j;            // j counts 0,1,2,... N-1
  70.   register int jp;            // Reverse counter:
  71.                     // N, N/2, N/4, N/4+N/2, N/8, N/8+N/2,
  72.                     // with N identical to 0 (mod N)
  73.   register short * tablj = index_conversion_table;
  74.   register int k;
  75.  
  76.   for(j=0, jp=0; *tablj++ = jp, ++j < N;)
  77.   {
  78.     for(k=N/2; k <= jp; jp -= k, k /= 2)
  79.       ;
  80.     jp += k;
  81.   }
  82. }
  83.  
  84.  
  85.                 // Fill in the W array
  86.                 // W[j] = exp( -I 2pi/N * j ), j=0..N-1
  87. void FFT::fill_in_W()
  88. {
  89.   register Complex * wj = W;
  90.   register int j;
  91.   for(j=0; j<N; j++)
  92.   {
  93.     register Complex arg = Complex(0,-2*PI/N * j);
  94.     *wj++ = exp(arg);
  95.   }
  96. }
  97.